/*
 *  SurfaceBox.cpp
 *  EPPI-0.0
 *
 *  Created by Roberto Ortega on 22/06/11
 *
 */

#include "surfacebox.h"

// ***
SurfaceBox::SurfaceBox()
:Surface()
{}

// ***
SurfaceBox::~SurfaceBox()
{}

// ***
SurfaceBox::SurfaceBox(ifstream &in, const string& na)
:Surface(na)
{
	dx = 1.0;
	dy = 1.0;
	dz = 1.0;
	
	nn = 1;
	mm = 1;
	ll = 1;
	
	string keyword;

	while (in >> keyword && keyword.compare("ENDSURFACE") )
	{		
		if (!keyword.compare("DIMS"))
		{
			in >> dx;
			in >> dy;
			in >> dz;
			
			cout << keyword << endl;
			cout << dx <<"\t" << dy << "\t" << dz <<"\t" << endl;
		}
		if (!keyword.compare("GRID"))
		{
			in >> nn;
			in >> mm;
			in >> ll;
			
			cout << keyword << "\t";
			cout << nn << "\t" << mm << "\t" << ll <<endl;
		}
		if (!keyword.compare("CONTACTPOINTS"))
		{
			in >> pointsOption;
			
			if (!pointsOption.compare("all"))
			{
				contactAllOption = true;
				
				cout << keyword << "\t";
				cout << pointsOption << endl;
			}
			
			if (!pointsOption.compare("file"))
			{
				in >> contactFileName;
				contactFileOption = true;
				
				cout << keyword << "\t";
				cout << pointsOption << "\t" << contactFileName << endl; 
			}
		}		
	}
	
	createBox();
	this->createContactPoints();
}

// ***
void SurfaceBox::createBox(void)
{
	int i;
	int j;
	
	double x;
	double y;
	double z;
	
	double u;
    double v;
	
	int elem;
    int node;
	
	vector<double> vCoordinates;
    vector<int> vElement;
	
	vCoordinates.resize(3);
	vElement.resize(4);
	
	// Plane z = +dz / 2
	for (j = 0; j < mm + 1; ++j)
	{
		for (i = 0; i < nn + 1; ++i)
		{
			u = dx / nn * i;
			v = dy / mm * j;
			
			//node = i + j * (nn + 1);
			
			x = -dx / 2 + u;
			y = -dy / 2 + v;
			z = +dz / 2;
			
			vCoordinates[0] = x;
			vCoordinates[1] = y;
			vCoordinates[2] = z;
			
			vLocalPoints.push_back(vCoordinates);
		}
	}

	for (j = 0; j < mm; ++j)
	{
		for (i = 0; i < nn; ++i)
		{
			elem = i + j * nn;
			
			vElement[0] = elem + j;
			vElement[1] = elem + j + 1;
			vElement[2] = elem + j + nn + 2;
			vElement[3] = elem + j + nn + 1;
			
			vLocalElements.push_back(vElement);
		}
	}
	node = vLocalPoints.size();
	
	// Plane z = -dz / 2
	for (j = 0; j < mm + 1; ++j)
	{
		for (i = 0; i < nn + 1; ++i)
		{
			u = dx / nn * i;
			v = dy / mm * j;
			
			//node = i + j * (nn + 1);
			
			x = -dx / 2 + u;
			y = -dy / 2 + v;
			z = -dz / 2;
			
			vCoordinates[0] = x;
			vCoordinates[1] = y;
			vCoordinates[2] = z;
			
			vLocalPoints.push_back(vCoordinates);
		}
	}
	
	for (j = 0; j < mm; ++j)
	{
		for (i = 0; i < nn; ++i)
		{
			elem = i + j * nn + node;
			
			vElement[0] = elem + j;
			vElement[1] = elem + j + 1;
			vElement[2] = elem + j + nn + 2;
			vElement[3] = elem + j + nn + 1;
			
			vLocalElements.push_back(vElement);
		}
	}
	node = vLocalPoints.size();

	// Plane y = +dy / 2
	for (j = 0; j < ll + 1; ++j)
	{
		for (i = 0; i < nn + 1; ++i)
		{
			u = dx / nn * i;
			v = dz / ll * j;
			
			//node = i + j * (nn + 1);
			
			x = -dx / 2 + u;
			y = +dy / 2;
			z = -dz / 2 + v;
			
			vCoordinates[0] = x;
			vCoordinates[1] = y;
			vCoordinates[2] = z;
			
			vLocalPoints.push_back(vCoordinates);
		}
	}
	
	for (j = 0; j < ll; ++j)
	{
		for (i = 0; i < nn; ++i)
		{
			elem = i + j * nn + node;
			
			vElement[0] = elem + j;
			vElement[1] = elem + j + 1;
			vElement[2] = elem + j + nn + 2;
			vElement[3] = elem + j + nn + 1;
			
			vLocalElements.push_back(vElement);
		}
	}	
	node = vLocalPoints.size();
	
	// Plane y = -dy / 2
	for (j = 0; j < ll + 1; ++j)
	{
		for (i = 0; i < nn + 1; ++i)
		{
			u = dx / nn * i;
			v = dz / ll * j;
			
			//node = i + j * (nn + 1);
			
			x = -dx / 2 + u;
			y = -dy / 2;
			z = -dz / 2 + v;
			
			vCoordinates[0] = x;
			vCoordinates[1] = y;
			vCoordinates[2] = z;
			
			vLocalPoints.push_back(vCoordinates);
		}
	}
	
	for (j = 0; j < ll; ++j)
	{
		for (i = 0; i < nn; ++i)
		{
			elem = i + j * nn + node;
			
			vElement[0] = elem + j;
			vElement[1] = elem + j + 1;
			vElement[2] = elem + j + nn + 2;
			vElement[3] = elem + j + nn + 1;
			
			vLocalElements.push_back(vElement);
		}
	}	
	node = vLocalPoints.size();
	
	// Plane x = +dx / 2
	for (j = 0; j < ll + 1; ++j)
	{
		for (i = 0; i < mm + 1; ++i)
		{
			u = dy / mm * i;
			v = dz / ll * j;
			
			//node = i + j * (mm + 1);
			
			x = +dx / 2;
			y = -dy / 2 + u;
			z = -dz / 2 + v;
			
			vCoordinates[0] = x;
			vCoordinates[1] = y;
			vCoordinates[2] = z;
			
			vLocalPoints.push_back(vCoordinates);
		}
	}
	
	for (j = 0; j < ll; ++j)
	{
		for (i = 0; i < mm; ++i)
		{
			elem = i + j * mm + node;
			
			vElement[0] = elem + j;
			vElement[1] = elem + j + 1;
			vElement[2] = elem + j + mm + 2;
			vElement[3] = elem + j + mm + 1;
			
			vLocalElements.push_back(vElement);
		}
	}	
	node = vLocalPoints.size();
	
	// Plane x = -dx / 2
	for (j = 0; j < ll + 1; ++j)
	{
		for (i = 0; i < mm + 1; ++i)
		{
			u = dy / mm * i;
			v = dz / ll * j;
			
			//node = i + j * (mm + 1);
			
			x = -dx / 2;
			y = -dy / 2 + u;
			z = -dz / 2 + v;
			
			vCoordinates[0] = x;
			vCoordinates[1] = y;
			vCoordinates[2] = z;
			
			vLocalPoints.push_back(vCoordinates);
		}
	}
	
	for (j = 0; j < ll; ++j)
	{
		for (i = 0; i < mm; ++i)
		{
			elem = i + j * mm + node;
			
			vElement[0] = elem + j;
			vElement[1] = elem + j + 1;
			vElement[2] = elem + j + mm + 2;
			vElement[3] = elem + j + mm + 1;
			
			vLocalElements.push_back(vElement);
		}
	}	
	
	// nLocalPoints and nLocalElements
	nLocalPoints = vLocalPoints.size();
	nLocalElements = vLocalElements.size();
}

// ***
double SurfaceBox::evalS(double& x, double& y, double& z)
{
	double aa = 0.05;
	double bb = 0.05;
	
	double rx = dx / 2;
	double ry = dy / 2;
	double rz = dz / 2;	

	S = pow(pow(fabs(x / rx), 2 / aa) + pow(fabs(y / ry), 2 / aa), aa / bb);
	S+= pow(fabs(z / rz), 2 / bb) - 1.0;
	
	return S;
}

//***
Eigen::RowVector3d SurfaceBox::evaldS(double& x, double& y, double& z)
{	
	double aa = 0.05;
	double bb = 0.05;

	double rx = dx / 2;
	double ry = dy / 2;
	double rz = dz / 2;
	
	dS.x() = pow(pow(fabs(x / rx), 2 / aa) + pow(fabs(y / ry), 2 / aa), aa / bb - 1);
	dS.x()*= sgn(x) * pow(fabs(x / rx), 2 / aa - 1) * 2 / (bb * rx);
	
	dS.y() = pow(pow(fabs(x / rx), 2 / aa) + pow(fabs(y / ry), 2 / aa), aa / bb - 1);
	dS.y()*= sgn(y) * pow(fabs(y / ry), 2 / aa - 1) * 2 / (bb * ry);
	
	dS.z() = sgn(z) * pow(fabs(z / rz), 2 / bb - 1) * 2 / (bb * rz);
	
	return dS;
}

//***
Eigen::Matrix3d SurfaceBox::evalddS(double& x, double& y, double& z)
{
	return ddS;
}